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Abstract 

Poikilothermic disease vectors can respond to altered climates through spatial changes in both 
population size and phenology. Quantitative descriptors to characterize, analyze and visualize 
these dynamic responses are lacking, particularly across large spatial domains. In order to 
demonstrate the value of a spatially explicit, dynamic modeling approach, we assessed spatial 
changes in the population dynamics of Ixodes scapularis, the Lyme disease vector, using a 
temperature -forced population model simulated across a grid of 4 x 4 km cells covering the 
eastern United States, using both modeled (Weather Research and Forecasting (WRF) 3.2.1) 
baseline/current (2001-2004) and projected (Representative Concentration Pathway (RCP) 4.5 
and RCP 8.5; 2057-2059) climate data. Ten dynamic population features (DPFs) were derived 
from simulated populations and analyzed spatially to characterize the regional population response 
to current and future climate across the domain. Each DPF under the current climate was assessed 
for its ability to discriminate observed Lyme disease risk and known vector presence/absence, 



© 2013 by the authors; hcensee MDPI, Basel, Switzerland. 

This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution license (http:// 
creativecommons.org/hcenses/by/3.0/). 

Author to whom correspondence should be addressed; justin.remais@emory.edu; Tel.: +1-404-712-8908; Fax: +1-404-727-8744. 
Conflict of Interest 

The authors declare no conflict of interest. 



Dhingra et al. 



Page 2 



using data from the US Centers for Disease Control and Prevention. Peak vector population and 
month of peak vector population were the DPFs that performed best as predictors of current Lyme 
disease risk. When examined under baseline and projected climate scenarios, the spatial and 
temporal distributions of DPFs shift and the seasonal cycle of key questing life stages is 
compressed under some scenarios. Our results demonstrate the utility of spatial characterization, 
analysis and visualization of dynamic population responses — including altered phenology — of 
disease vectors to altered climate. 
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1. Introduction 

Understanding the ecological response to anthropogenic environmental changes, including 
changes in climate, land use, land cover and other factors, requires quantitative tools to 
characterize, analyze and visualize dynamic changes in the populations of key organisms of 
interest. Developing such tools is made difficult by the fact that population responses to 
environmental change are spatially and temporally complex, particularly for organisms with 
multiple environmental life stages, such as those that participate in the transmission of 
vector-borne diseases (VBD). Disease vector populations may exhibit variations in seasonal 
timing and duration, and their generally non-linear response to environmental signals makes 
prediction of the risk posed by VBD under altered environmental conditions challenging [1- 
3]. When exposed to changing climatic conditions, vector distribution and the risks of VBD 
may shift substantially across time and space [1,4,5]. Yet, a great deal of uncertainty 
remains for many VBD systems [4,6-8], and little is known regarding the dynamic nature of 
the population response to climate change, particularly vector phenology (timing of life 
stages), seasonality and the duration of key population events. 

While some ecological analyses have characterized the dynamic population response of 
various plant and arthropod species to external forcings in a spatially explicit fashion (e.g., 
[9-12]), much analogous work on VBD has neglected the spatial domain [13,14]. Still, other 
work forgoes system dynamics, instead investigating the spatial patterns of static population 
measures, such as presence/absence or mean abundance (see, for instance, [15] for Lyme 
disease and [16] for hantavirus). Such analyses make use of statistical relationships between 
climate and habitat suitability to estimate, for instance, the potential changes in the 
distribution of habitat suitability for, or nymphal density of, Ixodes scapularis, the vector of 
Lyme disease [17,18]. This approach offers little insight into the nature of the population's 
response over time, such as shifts in peak population timing or variability in population 
density during key exposure periods (e.g., high season for recreational activities). Given the 
substantial and continuing disagreement regarding how chmate may change the distribution 
of VBD (e.g., [4,7]), analyses capable of assessing the relationship between exogenous 
forcings and population dynamics in space and time may provide such insights. 
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What is more, geovisualization of the dynamic VBD response to environmental change 
could provide key information (e.g., maps summarizing complex spatio-temporal 
phenomena) for developing policies to respond to shifting risk. Thus, geospatial tools for 
characterizing, analyzing and predicting the response of VBD to future changes are 
desirable, and these should emphasize dynamic phenomena known to be important for 
understanding risk, such as vector phenology and seasonality. Phenology — the timing of life 
stages — is known to be sensitive to climatic change and is an important determinant of the 
spatial distribution of arthropods [19,20]. Current models investigating arthropod 
distribution under future climates generally ignore phenology, instead, establishing a 
relationship between a vector's current abundance and key habitat characteristics and, then, 
applying that model to projected future conditions [21,22]. An examination of an organism's 
phenological response can reveal important, but subtle, impacts of changing climate. For 
instance, the date of flowering and fruiting have been shown to be important determinants of 
aspen distribution [23], and the date of first oviposition has been shown to be important for 
gypsy moth distribution [24]. Characterization of life stage-specific dynamic responses can 
highlight such subtle determinants of the distribution of vectors under the future climate. 

The seasonality of events may also shift under future conditions, with important 
consequences for VBD risk. For instance, vector populations may peak at certain times of 
the year, with peak incidence of disease occurring at other times (e.g., see [25] for Lyme 
disease). Some models of VBD response to climate change attempt to roughly characterize 
changes in seasonality (e.g., [26]); some integrate seasonal elements, such as temperature, 
humidity and daytime hours, through degree-day models (e.g., [27-30]), and still, others do 
not explicitly account for seasonality (e.g., [31]). A more detailed spatial representation of 
seasonal shifts would make it possible to characterize the potentially profound effect that 
environmental change may have on the length and timing of VBD transmission seasons. 

Here, we develop a spatially-explicit modeling approach for investigating the dynamic 
population responses of a disease vector of interest, with the goal of enhancing our 
understanding of future VBD risk. We introduce the concept of dynamic population features 
(DPFs), which provide information on population cycling, seasonal timing and phenological 
events across vector life stages. Importantly, we describe how analysis of such features — 
such as number and timing of population peaks (Table 1) — may be used to predict disease 
risk. 

To demonstrate the utility of this modeling approach, we examine the responses of the 
black-legged deer tick {Ixodes scapularis), the vector for Lyme disease, to changes in 
temperature across the eastern United States. /. scapularis is an excellent model organism 
with which to examine the influence of climate change on phenological and seasonal 
characteristics: it is known to be highly sensitive to environmental conditions, including 
temperature [6,15,25,32]. Furthermore, the three /. scapularis life stages (larva, nymph and 
adult) require different temperature conditions to support host finding or progress to the next 
life stage [28]. 

We explore the Lyme disease system as a case study, simulating /. scapularis population 
dynamics over the eastern United States using modeled climate data, and spatially 
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characterize, analyze and visualize key DPFs for each tick life stage. We examine DPFs 
from simulated dynamics under current climate conditions and compare these to observed 
data to ascertain which features best predict current levels of disease risk. We then project 
DPFs under two future climate scenarios and provide key geovisualizations of projected 
vector dynamics over the spatial range. We show, by characterizing and visualizing DPFs, 
how we can determine which population features best predict disease risk under current 
conditions and can then explore how future conditions may lead to shifts in these same DPFs 
in the future. We analyze DPFs in the context of /. scapularis and Lyme disease risk, but 
note that the approach shows promise for other organisms and disease systems. 

2. Methods 

2.1. Modeling Methodology 

The overarching analysis involved four key steps. First, a deterministic, dynamic population 
model was run, in parallel, over a large geographic area to generate spatially explicit 
simulations of population density in response to temperature variation. A daily time step was 
used in conjunction with the smallest grid cell size for which temperature data were 
available from a global circulation model. Second, simulated population dynamics were 
recorded at each grid cell for each vector life stage under current and future climate 
scenarios, and these were characterized in terms of their dynamic population features 
(DPFs), which were chosen to highlight population trends, seasonality or a combination of 
both. Third, DPFs were evaluated for their ability to predict the current distribution of 
vectors or human disease risk, using publically available data. Finally, DPF values found to 
be important determinants of current vector or disease distributions were visualized across 
the spatial domain for a range of future climate scenarios. We describe each of these steps in 
detail, with the application to Lyme disease, next. 

2.2. Lyme Model 

A twelve-stage temperature-driven Ufe cycle model of black-legged deer ticks (7. scapularis) 
(described in [28]) was adapted for high-performance simulation using Simulink (v. 7.0) and 
Matlab (R 201 lb) and executed across a cluster of 48 nodes. The model was run using daily 
(rather than monthly) mean temperatures, and the model was coded to be spatially explicit, 
executed in parallel at each grid cell across a large spatial domain. Temperature drivers 
shape simulated tick populations through degree-day functions, which model development 
delays, and through temperature-dependent activity parameters, which model host-seeking 
behaviors (see Supplementary Information, Table SI). Using a daily time step, model 
simulations at each cell were carried out over the domain under both the baseline and 
projected climate periods. Time series outputs, recorded daily, included questing adult (QA), 
questing nymph (QN), and questing larva (QL) populations. Given the rather large spatial 
resolution simulated from the perspective of Ixodes ecology (e.g., [34,35]), populations do 
not interact between grid cells {i.e., im/emigration are not modeled). 

2.3. Domain and Climate Data 

Simulation and subsequent analyses were conducted in the eastern United States across a 
domain of 4 x 4 km grid cells. Customized climate data across this grid were obtained from 
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the Regional Climate Model (Weather Research and Forecasting (WRF) 3.2.1) simulated at 
the University of Tennessee/Oak Ridge National Laboratory [36]. Daily temperatures at 
each cell were calculated as an average of daily minimum and maximum temperatures 
produced by climate simulations for a baseline time period (2001-2004) and two projected 
scenarios of differing severity (for 2057-2059). Projected scenarios. Representative 
Concentration Pathway (RCP) 4.5 and RCP 8.5, correspond to a continuous rise in radiative 
forcing to 4.5 W/m^ (moderate scenario) and 8.5 W/m^ (severe scenario), respectively, in 
2100 [37]. 

2.4. Dynamic Features of Ixodes Population Response to Seasonality Shifts 

Dynamic population features, chosen to highlight Ixodes population phenology and 
seasonality, were determined as described in Table 1 for each year. These were used to 
compare simulated population dynamics for each life stage at the grid cell level for three 
years of simulation under both the baseline and projected climate conditions. With the 
exceptions of the Mean or Median, calculated as the three-year mean or median of the 
simulated daily population, DPFs were calculated for each year of the simulation period and 
averaged to produce a final DPP value at each cell for each climate scenario. 

Population response DPFs included three-year Mean and three-year Median populations; the 
maximum population during each year (Peak Population); and the average number of local 
maxima per year (Peaks per Year). A 90-day moving window was used to identify each 
local peak through the course of the year. In order to identify peaks in the first 90 days of the 
simulation output, the last 89 days of simulation output were prepended to the output to 
provide a 90-day window. Similarly, the first 89 days of simulation output were appended to 
the simulation in order to aid identification of peaks in the final 90 days of the simulation. 

Seasonality DPFs, which were defined for each life stage, included two quantifications of 
season length and one of season timing. One classification of season length, termed IP to IP, 
was defined as the number of days between inflection points on either side of the annual 
maximum population. Inflection points were defined as changes in the concavity of the 
loess-smoothed (30-day window) population time series; changes in concavity were 
determined using a 3-point central difference equation on the smoothed time series. The 
second season length quantification was defined as the number of days from the annual 
maximum population to the annual minimum population, termed Peak to Trough. Wave 
Angle may be understood as the relative timing of each cell's season. To determine Wave 
Angle, continuous wavelet analysis was carried out using a Morlet wavelet [38] for the 
period of maximum power, -90.5 days (see Supplementary Information). DPFs combining 
seasonality and absolute population included IP Pop is the number of tick-days during each 
life stage's season (that is, the summation of the tick population for all days included in the 
IP to IP calculation). Additionally, UQ/IQR is used to estimate the period within the 
calendar year where the highest quartile simulated populations occur. Thus, UQ/IQR is 
defined by selecting the time points (days) in which the upper quartile populations occur, 
then taking the mean of the interquartile range of these time points (Table 1). 
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2.5. Comparison of DPFs to Observed Data 

DPFs obtained from the model as described above were fit to observed county-level /. 
scapularis presence (coded in three levels as absent, reported and established) and Lyme 
disease incidence (coded in four levels as none/minimal, low, medium and high) obtained 
from the Centers for Disease Control and Prevention (CDC) [33,39]. For all analyses, the 
four reported classifications were grouped into all possible dichotomizations (e.g., for Lyme 
disease, dichotomizations included minimal/none vs. low, medium and high; minimal/none 
and low vs. medium and high; and minimal/none, low and medium vs. high). DPFs were 
spatially averaged to the county level and compared to the observed (CDC) data using both 
area under the receiver operating characteristic curve (AUC) and logistic regression to 
ascertain each DPF's predictive ability. AUC (range: 0 to 1) quantizes model predictive 
accuracy for a dichotomous outcome, where a value of 0.5 indicates no predictive ability, a 
value of 1 indicates perfect discrimination and a value of 0 indicates lack of discrimination. 
To assess potential spatial variation in the ability of DPFs to predict Lyme disease risk, 
AUCs for selected DPFs were also determined for counties in three regions (Northeast, 
South and Midwest), as defined by the US census. 

2.6. Spatial Sensitivity Analysis 

To assess DPF sensitivity to spatial autocorrelation, the above-defined model outputs were 
also fit to observed Lyme disease and tick prevalence data using both a (non-spatial) logistic 
model and a spatial logistic model. The spatial regression model is defined as follows: 

logit[P{yi=l)\=a+\i+3Xi 

where y,- is the dichotomized observed Lyme disease category, a is the overall baseline risk, 
\i is the county-specific spatial random effects and (3 represents the log odds ratio associated 
with DPFs of population response (X,). The model used controls for the effects of spatial 
autocorrelation using an intrinsic conditional autoregressive (CAR) model [40]. The CAR 
model, often formulated by the conditional distribution of "kj, given its neighbors, assumes 
that Xi, for each county, ;, is a spatial average of its neighbors. The conditional distribution is 

Gaussian, with mean ZT^,^,/"^' and variance x^/m,, where denotes that county, 

Till 

shares a boundary with county and m/ is the number of boundary-sharing neighbors for 
county /. The percent change in (3 and the magnitude of parameter, t^, which controls the 
degree of spatial similarly, were used to assess the improvement in fit provided by the 
addition of the spatial term to the logistic regression model. 

3. Results 

Simulated population dynamics for the three questing tick life stages were obtained from a 
temperature-forced model simulated across a grid of 4 x 4 km cells covering the eastern 
United States under both baseline and projected climate. Based on these simulated 
population dynamics, ten dynamic population features (DPFs) were derived and analyzed 
spatially to characterize the population response to current and future climate across the 
domain. Pairwise correlations between DPFs were determined, and each DPF under the 
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current climate was assessed for its ability to discriminate Lyme disease risk and vector 
presence/absence using observed data from the US CDC. 

3.1. Correlation among DPFs 

To determine the degree to which DPFs are collinear and, thus, to quantify related aspects of 
the population curve, correlations were calculated for all pairs of DPFs. There was strong 
correlation between two different DPFs of absolute population {Mean and Median) in all 
questing life stages over the spatial extent of analysis (rg > 0.96; see Table 2). Peaks per 
Year and IP Pop were also highly correlated (r^ > 0.88) with absolute population DPFs and 
each other in the questing adult stage, while Peaks per Year was also highly correlated {r^ > 
0.90), with absolute population DPFs and IP to IP in the questing nymph and larval stages. 
Only Wave Angle was inversely correlated with all other DPFs in the nymph and adult 
stages. With the exception of its correlation with UQ/IQR (QA: rg = -0.77; QN: = -0.81), 
this inverse relationship was weak (rj > -0.63). In the questing larval stage, only IP to IP 
showed a weak inverse correlation with other DPFs (save for its correlation with UQ/IQR). 
No timing DPF was strongly associated with any other timing DPF in any life stage. 

3.2. Comparison of DPFs to Observed Data 

Each DPF was fit to observed tick presence or to Lyme disease risk from the US CDC, 
producing an AUC value, a discriminatory index that allows comparison of continuous 
predictions to dichotomous observations without requiring subjective cut points. Peak 
Month and Peak Population showed the greatest discriminatory ability across all life stages 
when compared with CDC Lyme disease risk data (Table 3; AUC = 0.54 to 0.90). Among 
most DPFs, the minimal vs. high or the minimal/low/medium vs. high dichotomizations gave 
higher AUC values than the remaining dichotomizations. The dichotomization of minimal 
vs. low/medium/high performed worst in AUC analyses across all questing life stages. Peak 
to Trough, IP to IP, IP Pop, UQ/IQR, and Wave Angle showed inconsistent predictive 
ability over all three life stages. Of these, only IP Pop in the QN life stage showed some 
improved discriminatory ability. 

Comparison of DPFs to CDC tick presence data showed markedly less discriminatory 
ability. Peak Population, Peak Month, Peak to Trough and UQ/IQR were statistically 
significant predictors of tick presence across all questing life stages and dichotomizations of 
CDC data, but AUC values were low and ranged from 0.53 to 0.71 for QA, from 0.54 to 
0.69 for QN, and from 0.52 to 0.69 for QL. For QL Wave Angle and tick presence, AUC 
ranged from 0.66 to 0.7. Along with QA Peak Month and QA Peak Population, this was 
among the strongest predictors of tick presence. However, for the dichotomization minimal/ 
none vs. low/medium/high, DPFs, Mean, Median, Peaks per Year, IP to IP, IP Pop and UQ/ 
IQR, were uniformly non-significant across all questing life stages. Peak Month AUC values 
were very similar for QN and QL, and slightly lower for QA. 

3.2.1 . Regional Analyses — A regional analysis, analogous to the preceding AUC 
analysis, was carried out for Peak Month and Peak Population, to ascertain the degree to 
which discriminatory ability varies by location. In most cases, AUCs for Peak Month and 
Peak Population were observed to be highest (Table 4) in the Midwest region of the eastern 
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United States (AUC: 0.80 to 0.96), where the AUC was most often statistically significant. 
AUC values in the North were consistently lower than the Midwest (AUC: 0.71 to 0.78), for 
both DPFs. In the South, Peak Month usually demonstrated higher statistically significant 
predictive ability for QN and QL than Peak Population, while both DPFs demonstrated 
similar predictive ability for QA. 

3.2.2. Spatial Sensitivity Analysis — Regression coefficients produced in conditional 
autoregressive models for each DPP at each life stage on CDC data were not substantively 
different from those produced by a non-spatial regression (not shown). Conditional 
autoregressive models exhibited low values of (0.003, 0.06) relative to the intercept and/or 
the parameter value, indicating that the contribution of spatial autocorrelation is small. 

3.3. Shifts in Geographic Distribution of DPFs 

DPFs in all cells (N = 262,152) under future climate scenarios were significantly (p < 0.001) 
different than those from baseline simulations. Simulated mean, median and peak 
populations all show increases across most of the eastern United States, with the largest 
increases in the RCP 8.5 scenario. In particular, 18.8%, 7.7% and 4.1% of cells showed an 
increase of an order of magnitude or greater in peak QA, QN and QL population, 
respectively, in the RCP 8.5 scenario. Relative to baseline, regions of highest projected 
mean, median and peak tick population expanded both northward and southward to encroach 
upon the areas of low DPP values occurring across the Appalachian mountain range (e.g.. 
Figure 1). 

While both projected scenarios showed simulated questing life stage average Peaks per Year 
that were significantly different from the baseline case and the two projected scenarios were 
significantly different from one another, there was no substantive change (<0.01 peaks per 
year) when these comparisons were made across the entire domain. Variations in population 
response between scenarios, as shown by Peaks per Year, demonstrate the lack of uniform 
response (Figure 1). The number of cells projected to have two or more population peaks per 
year for the QA life stage increased in northern regions, while there was a net decrease in 
the number of cells with more than one peak in southern regions. In contrast, for the QL and 
QN life stages, the number of cells experiencing more than one population peak per year 
decreased uniformly across the spatial domain. 

Under scenarios RCP 4.5 and RCP 8.5 compared to baseline, the season length for QA, as 
defined by IP to IP, has a pronounced increase from 50-70 day seasons to 60-90 and 
80-100 day seasons, respectively, in the South. In the North, there is a net decrease in QA 
season length for both future scenarios (Figure 1). For QN, projected season length remained 
the same as baseline, though there was a slight decrease in season length for much of the 
northern portion of the study area in both projected scenarios, and an increase of 
approximately 40 days in season length in a small portion of the land-locked Midwest. 
Projected QL season length showed a decrease of approximately 10 days in southern areas, 
while the overall geographic area with higher season lengths decreased with increasing 
scenario severity. Changes in projected season length, as defined by Peak to Trough, were 
approximately uniform across the domain (not shown). QA Peak to Trough lengthened on 
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average by eight days in RCP 4.5, and by 1.4 days in RCP 8.5, as compared to baseline. QN 
and QL Peak to Trough shortened by 56.6 and 31.1 days, respectively, in the RCP 8.5 
scenario, and by 15 and 10 days, respectively, in the RCP 4.5 scenario. 

The exposure DPF, IP Pop, which counts the number of tick-days during a season bounded 
by inflection points around the yearly maximum population, bears little similarity to IP to IP 
maps (Figure 1). IP Pop increases in exposure-time across the domain as the severity of the 
projected scenarios increases. The regions with the highest baseline number of ticks, present 
during the inflection point defined season, spread outward in both the north and south 
direction and center around the Midwest and the Northeast. 

Wave Angle results (see Supplementary Information; Figure SI) showed that the dynamics 
of all life stages under both projected climate conditions lag behind those at baseline climate 
by <4 days across the simulated domain. The projected month of peak population (Peak 
Month) and UQ/IQR for QA generally shifted to earlier months. However, with increasing 
scenario severity, QN and QL generally shifted to later months across the geographic area. 

Wave angle results (see Figure SI) showed that the dynamics of all life stages under both 
projected climate conditions lag behind those at baseline climate by <4 days across the 
simulated domain. The projected month of peak population (Peak Month) and UQ/IQR for 
QA generally shifted to earlier months. However, with increasing scenario severity, QN and 
QL generally shifted to later months across the geographic area. 

4. Discussion and Conclusions 

When examining the response of vector populations to climate change, shifts in phenology, 
seasonality and other dynamic characteristics can be anticipated across the spatial range and 
life stages of the organism of interest. Risk of VBD is dependent on both timing and 
probability of exposure to the vector, and thus, characterizing the dynamic population 
response over space is crucial in order to anticipate and manage potential future risks. Here, 
we provide a framework for evaluating both static and dynamic effects of climate change on 
populations over large geographic areas, using spatially explicit simulation of a climate- 
driven, stage-structured population model. 

Our findings with respect to Ixodes scapularis illustrate both the methodology and its utility. 
The derivation and analysis of dynamic population features are key to the analytical 
approach. DPFs provide quantitative information about a range of population characteristics 
and allow for comparison between dynamic simulation output and observed disease data, as 
well as between baseline and projected climates. Absolute DPFs, such as Mean, Median and 
Peak Population can be interpreted as indicators of survivorship, while timing DPFs, such as 
number of days from the yearly maximum population to the yearly minimum population 
{Peak to Trough) and month in which peak population occurs (Peak Month) characterize the 
timing and length of a given life stage's season. 

In the case of Ixodes, DPFs associated with the peak of the simulated population curve. Peak 
Population and Peak Month, proved to be the most important in predicting high risk of 
Lyme disease, though all DPFs showed some level of discriminatory ability. AUC analyses 
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showed that dichotomizations isolating high risk improved discriminatory ability across all 
DPFs and life stages. Aggregation of medium and high risk also showed improved 
discriminatory ability across life stages and DPFs as compared to the minimal vs. low/ 
medium/high dichotomization. This trend of improvement, as high disease risk is 
progressively isolated into a single category, suggests strongly that these DPFs are useful in 
predicting the timing and location of higher Lyme disease risks. 

When DPFs are examined for two projected climate scenarios, we show that the dynamic 
population response of /. scapularis is not uniform across life stages and varies over space. 
Spatial shifts in temporal features include geographic shifts in season, and these shifts are 
not consistently northward as one might intuitively hypothesize. While the month in which 
the greatest number of ticks are questing (Peak Month) is delayed for the adult life stage 
(Figure 1), QN and QL peaks do not show geographically uniform shifts to earlier questing 
season. Also, Peak to Trough and IP to IP indicate potential changes in season length in 
projected scenarios. Spatial shifts in absolute DPFs, such as Peak Population, vary by 
region. For instance, the peak populations in the Midwest and the Northeast regions are both 
expected to rise far more as compared to the Appalachian mountain range or the Gulf Coast, 
where these populations are expected to remain more stable. 

Although the finding that QL Peak Month and Peak Population show high predictive ability 
for Lyme disease risk is significant, the causal implications of this finding, and others like it, 
must be interpreted cautiously. Disease risk is not directly related to the questing larval 
stage, which takes the first blood meal in the lifecycle, and thus, is responsible for Lyme 
transmission only under the rare circumstance that larvae are infected transovarially. 
Likewise, QN Peak Month and Peak Population have similar AUC values for all 
dichotomizations of Lyme disease risk, an effect driven largely by the similarity of tick 
response to temperature in these two life stages, rather than mutual causal relationships with 
disease. Complex temporal relationships are inherent in these populations: questing nymphs 
and questing larvae, for instance, peak at approximately the same time of year, and their 
populations in a given location are ostensibly correlated, though the QN population does not 
result from the QL population in the same year, but rather previous years' QL. 

As in other ecological modeling analyses, data quality determines the utility of this analysis 
framework for a given system. In our analysis, CDC data quality may account for the lack of 
significant AUCs of DPFs in comparison to the observed tick data. Tick presence/absence 
data are collected using a variety of methods, such as dragging and deer surveys, often under 
serious resource constraints [33]. Rather than providing consistent, systematic information 
about tick presence and absence, the national tick dataset offers a coarse categorization 
derived from disparate information. This is in contrast to the national Lyme disease dataset, 
which is based on a consistent reporting standard. Given the higher quality of data collected, 
this dataset is more useful in substantiating the results of our model. 

Other climate factors besides temperature, such as humidity, have been shown to affect 
Ixodes spp. activity [41,42] and correlate with human Lyme disease risk [43]. The 
population model used here did not incorporate Ixodes' response to humidity, and although 
our simulated population data demonstrated good correspondence with Lyme incidence, it is 
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possible that including other key environmental variables may yield yet greater 
correspondence. Likewise, host and pathogen populations were not considered in our 
analysis, which was limited to vector dynamics. Relatively little research has been done on 
the potential population responses of Borrelia spp. under altered climate conditions. 
However, it has been suggested that changes in Ixodes phenology in response to climatic 
changes may affect the evolution of various tick-borne pathogens, so as to modify their 
lifespan, transmission and pathogenicity [44]. Host dynamics can also greatly impact 
infected vector density and consequent human risk in a variety of VBD systems [45^7]. In 
the case of Lyme disease, the abundance of key hosts, such as mice and chipmunks, has 
been shown to predict the density of infected nymphs in eastern deciduous forests [48]; in 
other areas, such as the southern United States, lizards are believed to exert a dampening 
effect on the spread of Lyme disease, due to poor host competence or zooprophylactic 
effects [49,50]. Including host, vector and pathogen dynamics in a combined model would 
pose significant methodological and computational challenges, but is also likely to add 
greatly to our mechanistic understanding of shifting VBD risk under future environmental 
conditions. We note that a similar simulation, summarization (e.g., DPFs) and analysis 
approach can be pursued with such a combined model; yet, other summarizations (e.g., Rq) 
become available for geovisualization in that context (e.g., [51,52]). 

The methodological contributions made by the modeling analysis described here are 
considerable. We provide a quantitative assessment of population dynamics — with potential 
consequences for disease risk — under future climates, which is made possible by use of a 
spatially-explicit, mechanistic model [53]. Our spatial characterization of DPFs allows a 
detailed visual assessment (e.g.. Figure 1), alongside a quantitative analysis, of the dynamic 
population response to future climate, revealing potential changes that are non-intuitive. For 
instance, across the eastern United States, under projected temperatures as compared to the 
baseline scenario, nymphs and larvae are projected to arrive at their peak population earlier 
in the season, while adults are projected to reach peak population later in the season (Figure 
1). The approach taken here also highlights the value of modeUng abundance, which, unlike 
habitat suitability or other static measures, allows for the examination of phenology and 
seasonality among life stages and the potential implications for (and correlation with) 
disease risk. For instance, IP to IP indicates that the length of larval "season" is stable across 
the three temperature scenarios, while the adult, and, to a lesser extent, the nymphal stages 
exhibit "seasons" that are strongly sensitive to the projected increasing temperatures. Such 
life-stage-specific responses in time and space would be unapparent using traditional 
methods that examine, for instance, aggregate, annual effects. 

We caution above against a causal interpretation of a DPF's predictive power. A strong 
correlation between a DPF and observed disease incidence may not represent a causal 
relationship, but such a finding can raise hypotheses that ultimately lead to greater 
mechanistic understanding of the relationship between vector populations and disease risk in 
space and time and, thus, an improved causal understanding. Finally, population models, 
such as the one examined here, can also be used to evaluate the efficacy and economy of 
potential public health interventions [53], such as vector or host control (e.g., [54-56] for 
Lyme disease). A coupled analysis of the effect of temperature in the presence of a vector 
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control program would be an obvious extension of the approach, and such an application of 
this model is possible for many different vectors, interventions and diseases. 

We have demonstrated the ability of a spatially-explicit dynamic population model to 
discriminate between dynamic population features most strongly associated with disease 
risk, as well as to characterize the geographically varied response of /. scapularis life stages 
to climate dynamics. Use of such an approach to describe shifts in dynamics is not limited to 
Lyme disease. The technique may provide new insights into the dynamic responses of a 
range of disease vectors to environmental changes, particularly shifts in their seasonal and 
phenological features. Such analyses may provide helpful information about the consequent 
risk of vector-borne disease under future conditions. 
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Figure 1. 

(A) Log of Peak Population, (B) IP to IP, (C) log of IP Pop, (D) Peak Month, (E) UQ/IQR 
and (F) Peak to Trough for questing adults (QA), questing nymphs (QN) and questing larvae 
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Table 1 

Dynamic population features (DPFs) of population response. 
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